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Abstract. The relaxation of the specific heat and the entropy to their equilibrium values is investigated 
numerically for the three-dimensional Coulomb glass at very low temperatures. The long time relaxation 
follows a stretched exponential function, fit) = fo exp [— (t/r) 1 , with the exponent (5 increasing with 
the temperature. The relaxation time follows an Arrhenius behavior divergence when T — > 0. A relation 
between the specific heat and the entropy in the long time regime is found. 



PACS. 61.43.Fs 
transition 



Glasses - 61.43.Dq Amorphous semiconductors, metals, and alloys - 64.70.Pf Glass 



> 
SO 
OS 
CO 

so 
O 
O 
O 

c3 



T3 

C 

O 
O 



x 



1 Introduction 

The relaxation in most glassy systems such as structural 
glasses, ionic conductor supercooled liquids, polymer, col- 
loid, and spin glasses |jj deviates strongly from a single 
exponential relaxation at some temperature T* well above 
the static transition. In the three dimensional Ising spin 
glass model Ogielski Q observed that the long time regime 
of the relaxation functions is well approximated by the 
non-exponential function 



f(t) = f t-*exp [-(t/rf] 



(1) 



below T* — Tq, where T G is the Griffiths temperature ||. 
In supercooled liquids || the characteristic features of the 
relaxation processes in the long time regime is captured 
by a stretched exponential or Kohlrausch- Williams- Watts 
decay function 



/(t) = / exp [-{t/rf 



(2) 



with < j3 < 1. This behavior can be obtained by a 
superposition of purely exponential relaxation process but 
subject to a broad distribution of relaxation times JfJ. 

The Coulomb glass is a prominent example of 

glassy systems. In heavily doped crystalline semiconduc- 
tors, amorphous semiconductor-metal alloys, and granular 
metals, it plays an important role as a semiclassical model 
for systems of localized states. The out of equilibrium dy- 
namical behavior of the Coulomb glass has been studied 
previously by Schreiber et al. H and Perez-Garrido et al. 
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[ |10[ . They determined numerically the transition probabil- 
ities between low-energy many-particle states, and stud- 
ied the eigenvalues of the transition probability matrix. 
A broad distribution of relaxation times over several or- 
ders of magnitude was found in both cases, which reflects 
the glassy behavior of this system. This broad distribution 
bring to a power law relaxation function |ro| . Moreover, 
Wappler et al. jll) used the damage-spreading algorithm 
to study the temporal evolution of the system in thermal 
equilibrium, and found evidence for a dynamical phase 
transition at a temperature T* which depended of the de- 
gree of disorder considered. Yu studied the time devel- 
opment of the Coulomb gap considering a self-consistent 
equation for the density of states [|l2| . She also observed 
that very long time scales are involved. Diaz-Sanchez et 
al. |T^,|l4]|, considering the dynamic in the configuration 
space, studied non-ergodic effects on the specific heat and 
found very long relaxation times. Pastor et al. pre- 
dicted the existence of a glass phase for a mean-field model 
of interacting spinless fermions in the presence of disorder. 
Experimental measurement in electronic system of an An- 
derson insulator jl6) showed a glass phase with aging phe- 
nomenon appearing. There it was found that the common 
relaxation law can be fitted by a stretched exponential 
function. 

The aim of this paper is to study the shape of the re- 
laxation process of the specific heat and the entropy in the 
Coulomb glass at very low temperatures and in thermal 
equilibrium. The paper is organized as follows: Section 2 
introduces the Coulomb glass model. Section 3 describes 
the numerical procedure, i.e. how the low-energy many- 
particle states are obtained numerically, and the dynamic 
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used in the simulation. In Section 4, we use this procedure 
for the investigation of the relaxation process of the spe- 
cific heat and the entropy. Finally, in Section 5 we extract 
some conclusions. 



from one up to four sites. In the second step this set S is 
improved by means of the thermal cycling method, which 
combines the Metropolis and local search algorithms. The 
third step completes the set S by systematically investi- 
gating the surroundings of the states previously found. 



2 The model 

A practical model to represent the Coulomb glass prob- 
lem for simulating an impurity band of localized states 
in lightly doped semiconductors, when quantum interfer- 
ence can be neglected, was proposed in references @J§- 
Later it was also applied to simulate granular metals ||17| 
and conducting polymers |Q . Following reference Jl!| , we 
consider a three dimensional cubic lattice with N sites. 
The sites can be occupied by KN (0 < K < 1) elec- 
trons which are interacting via an unscreened Coulomb 
potential. Background charges —K are added at each of 
the lattice sites, guaranteeing electro-neutrality. The dis- 
order is simulated by a random potential ej. Their values 
are uniformly distributed between —B/2 and B/2. This 
model is represented by the Hamiltonian 



E 



(ni-K)( nj -K) 



(3) 



where rii € {0, 1} denotes the occupation number of site 
i, is the distance between sites i and j according to 
periodic boundary conditions pb]| , and [i is the chemical 
potential. The lattice spacing is taken as unit of distance. 

In this paper, we focus in the half-filled impurity band 
K = 1/2, where there is particle-hole symmetry and the 
chemical potential /i equals zero at any temperature in 
the thermodynamic limit. We take B — 2 in this paper, 
similar results to the ones presented here are found for 
other values of B [B3. 



3.2 Dynamic in the configuration space 

We now present a method to study the influence of the 
duration of the measurement, i.e. the observation time 
r m , on the expectation value of the specific heat c and 
the entropy S Ji4| . This method defines a dynamic in the 
configuration space. 

During T m the state of the sample travels randomly 
through its configuration space. A transition between con- 
figurations / and J is allowed if the transition time tjj 
is shorter than r m . So, at finite r m only transitions with 
characteristic time tjj shorter than r m contribute to the 
specific heat. Thus, for a given r m , we consider two con- 
figurations as connected if their tjj is shorter than r m , 
and we group the configurations in clusters according to 
these connections. These clusters correspond to regions of 
the configurations space being isolated from each other on 
the time scale r m . 

The characteristic transition time between configura- 
tions I and J is [@,[IJ] , 



TIJ 



T Q exp (2 ra /a) exp (Eu/T) /Z . (4) 



In this equation, the quantity To is a constant of the or- 
der of the inverse of the phonon frequency, To ~ 10~ 13 s. 
The sum is the minimized sum over all hopping distances 
between sites which change their occupation in the tran- 
sition I —> J. a denotes the localization radius, Ejj = 
max(Ei, Ej) where Ej is the energy of the state /, T is 
the temperature, and Z is the partition function. 

Provided, at the beginning of the measuring process, 
the sample is in one of the states of the cluster a, we 
Our aim is to simulate the temporal evolution of the Coulombmeasure as value of the mean energy in this cluster, 
glass at very low temperatures in thermal equilibrium. For 



3 Numerical procedure 



that we first obtain a set S of low-energy many-particle 
configurations and later define a dynamic between these 
configurations. We apply this procedure to study the re- 
laxation of the specific heat and the entropy to their equi- 
librium values. 



3.1 Low-energy configurations 

We find the low-energy many-particle configurations by 
means of a three-steps algorithm Jl4| . This comprises lo- 
cal search 1 21 , 22]|, thermal cycling [|23[ , and construction 
of "neighboring" states by local rearrangements of the 
charges pl| , p2| . The efficiency of this algorithm is illus- 
trated in reference [ 14j . In the first step we create an ini- 
tial set S of metastable states. We start from states chosen 
at random and relax these states by a local search algo- 
rithm which ensures stability with respect to excitations 



(E) a = J2 E i p i . 



with 



Pi = exp 



-Ej 
T 



(5) 



(G) 



where Z a denotes the value of the the partition function 
of this cluster. The entropy of the cluster is given by 



lea 



Pi In (Pi) 



(7) 



Assuming that the clusters are in thermal equilibrium 
the specific heat for a cluster a of states is: 



(E 2 )a - (E)t 

T 2 N 



(8) 
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The expectation value of c, (c) , that is the mean value 
of repeated measurements, is given by the weighted aver- 
age of c a : 

(c(T,r m )> =]T Cct P Q , (9) 

a 

where P a is the probability to find the sample in one of the 
states belonging to the cluster a. Similarly the expectation 
entropy, (S), is given by 



(S(T,r m )> =^TS Q P Q 



(10) 



Thus (c) and (S) depend via the cluster structure and 
P a on r m , and also on T. Both the cluster structure and P a 
are influenced by the history of the sample. In our numer- 
ical experiments, we presume that the sample has reached 
thermal equilibrium before the measurements; thus P a = 
Z a /Z. 



4 Relaxation functions 

We study the process of relaxation of the specific heat 
and the entropy according to the methods presented in 
the previous section. To make the influence of r m directly 
visible, we consider the rate of the values of the specific 
heat and the entropy for finite and infinite (equilibrium) 
duration of measurement, respectively: 



q c (T,T m ) 



(c(r,r m )) _ (c(T,r m )) 
(c(T, oo)) c eq 

(S(T,r m )) _ (S(T,r m )) 
S (-■ , 



(11) 



(12) 



(5(T,oo)) u C q 

where c eq and S eq are the equilibrium values. Moreover, in 
order to characterize the shape of the relaxation process 
in the long time regime we study the following functions, 



R C (T, r m ) = 1 - q c (T, r m ) 
R S (T, r m ) = 1 - q B (T, r m ) 



(13) 
(14) 



which go to zero at the equilibrium. 

In our simulations wc consider the localization radius 
a equals always 0.2 |p4j . Investigating the physical proper- 
ties of macroscopic systems, we have calculated ensemble 
averages, and have compared the results for different sam- 
ple sizes. In order to ensure the convergence of the results 
on the number of low-energy many-particle configurations 
in the set S, we have taken from 25000 to 75000 configu- 
ration into account for the cases studied here (the number 
of configurations to consider depends on the temperature 
and system size). This election ensures that the width of 
the related energy interval exceeds the temperature by at 
least a factor of 25. As in reference for N > 512 the 
result are free of finite-size effects. Here we present the 
simulations for N — 512 taking into account 1000 samples 
for ensemble averages. In order to simplify the notation 
we consider the variable t — T m /ro, i.e. our unit of time is 

TO- 



Wc first study the shape of the relaxation functions R c 
and R s . Figure 1 shows R c and R s as a function of t, for 
five different temperatures. Here we can see that the long 
time regime, R c < 0.5 and R s < 0.5, may be reasonably 
fitted by a stretched exponential function. In the worse 
case, T — 0.014, the 75% of the stretched exponential 
function is compared with the data although for T — 0.006 
it becomes the 80%. For all temperatures considered here 
R c is equal to R s in the long time regime (within the 
errors). This form of the relaxation functions brings us to 
think in a broad distribution of relaxation times in the 
system, as it has been found in references p|,^0[. We have 
also attempted a fit of our results with equation (1) but 
the better fit is made with x ~ 0, i.e. we recover equation 
(2). 

The influence of the temperature on the relaxation pro- 
cess is studied. The values of f3, obtained by fitting R c and 
R s with equation (2), are presented in Figure 2. We see 
that increases with the temperature and it is the same 
for both relaxation functions, R c and R s , at each temper- 
ature (within the errors). This behavior of /? is also found 
in other glassy systems pEfl . We could expect some tem- 
perature T* , higher than the temperatures studied here, 
where ft = 1, recovering the exponential relaxation pro- 
cess at T > T* . Unfortunately our method does not let 
to investigate higher temperatures because we must take 
a finite number of configurations into account in the set 
S (we have taken the maximum of configurations that we 
have been able in order to make the calculations in a rea- 
sonable time). 

For an estimate of the long time relaxation we use r 
from the fit of equation (2) of R c and R s . As we can see in 
Figure 3, r r is reasonably fitted by an Arrhenius behavior 
divergence as T — ► 



T r (T) = aexp[T /T] 



(15) 



with a = 1.4 x 10 8 and T = 0.1106 for R c ; for R s we 
have a = 1.2 x 10 9 and T = 0.0946. The differences in a 
and To from both relaxation functions are because of the 
sensitivity of r on the details of the fit. 

We now consider the relation between the specific heat 
and the entropy in the long time regime. At the equilib- 
rium both are connected by the relation 



d(S) 
dT 



N I \ 



(16) 



As we will see below due to the temperature dependence 
of R c and R s this thermodynamic relation only hold at 
the equilibrium, i.e. when t 3> r r . We have seen that R c ks 
R s = R in the long time regime, so that q c rs q s — q and 
we can write from equations (11) and (12) 



(c) = q(T,t)c, 



-cq ? 



(17) 



and 



(S)=q(T,t)S eq , (18) 

after we have substituted (c) and {S) in equation (16) and 
regrouped terms we have 



fdS, 



cq 



V dT 



N \ dq_ 

rp C cq I q — ^cq qj, 



(19) 
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The left term of this equation is zero while the right term is 
different of zero when the system is not at the equilibrium. 
So that we can write a new relation for the entropy and 
the specific heat. It follows as 

9{S) dq N 

At the equilibrium dq/ dT vanishes and we recover equa- 
tion (16). In order to check this expression in Figure 4 we 
compare the specific heat calculated from this equation, 
and from the values of the entropy, with the specific heat 
calculated from equation (9). The derivatives are approx- 
imated with finite differences. We can see that both agree 
very well in the long time regime. 



5 Conclusion 

We have studied the relaxation process of the specific heat 
and the entropy for the three-dimensional Coulomb glass 
at very low temperatures and in thermal equilibrium. The 
long time relaxation regime follows a stretched exponen- 
tial function, with the exponent /3 < 1 increasing with the 
temperature, which is an indication of a broad distribu- 
tion of relaxation times in the system. From these results, 
we could expect a dynamical transition at some temper- 
ature T* (above of the temperatures studied here) where 
(3—1. The exponential relaxation process would be recov- 
ered for T > T*. In reference Q it was found T* > 0.03 
for the two-dimensional Coulomb glass although the ef- 
fects of distance-dependent transition probabilities were 
not taken into account. Moreover, the relaxation functions 
found here in the thermal equilibrium are very different 
from the power law relaxations found when the system 
is out of equilibrium flo|| . Nevertheless, in out of equilib- 
rium experiments stretched exponential relaxation func- 
tions are found jL6). The relaxation time diverges as an 
Arrhenius law when T — > 0. Finally, we have found a rela- 
tion between the specific heat and the entropy in the long 
time regime. 
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Fig. 3. Relaxation time r r as a function of the inverse of the 
Fig. 1. Relaxation functions of the specific heat as a function temperature, obtained by fitting R c (•) and R a (*) with the 
of t, for the temperatures T = 0.006 (o), 0.008 (■), 0.01 (o), function f(t) = f cxp [-(f/r r ) /3 ]. The lines are fits with the 
0.012(a) and 0.014 (*). Inset the relaxation functions of the function r r (T) = aexp[To/T]. The values of r r correspond to 
entropy for the same temperatures. The solid lines correspond the values obtained from the fit in Figure 1. 
to fits with the function f(t) — f exp [— (t/r)' 3 ] . 
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Fig. 2. Parameter /3 as a function of the temperature, ob- 
tained by fitting R c (•) and R 3 (o) with the function f(t) — 
/ exp[-(t/r)^]. 



Fig. 4. Specific heat as a function of time obtained from equa- 
tion (20) and the values of the entropy (*) and from equation 
(9) (•), for T = 0.008. 



